Set-up, Load Libraries

knitr::opts_chunk$set(echo = TRUE, warning = FALSE, eval = TRUE)
knitr::opts_knit$set(root.dir = normalizePath("../"))
#always looking from the perspective of where the rmd is, "../" says look outside of this folder
#~ is a shortcut for your home directory/root folder

#install.packages("OpenStreetMap")
library (OpenStreetMap) 
## Warning: package 'OpenStreetMap' was built under R version 3.3.3
library (rgdal)
## Loading required package: sp
## rgdal: version: 1.2-5, (SVN revision 648)
##  Geospatial Data Abstraction Library extensions to R successfully loaded
##  Loaded GDAL runtime: GDAL 2.0.1, released 2015/09/15
##  Path to GDAL shared files: C:/Users/ebola/Documents/R/win-library/3.3/rgdal/gdal
##  Loaded PROJ.4 runtime: Rel. 4.9.2, 08 September 2015, [PJ_VERSION: 492]
##  Path to PROJ.4 shared files: C:/Users/ebola/Documents/R/win-library/3.3/rgdal/proj
##  Linking to sp version: 1.2-4
library(sp) 
library(maptools)
## Checking rgeos availability: TRUE
library (raster)

Lab 4 Introduction to spatial data in R

In this lab you will begin to get used to using R with spatial data and practice your map making in R for a variety of data types.

Start by reading through sections 7.1-7.3 in Arnold and Tilton.2015. Humanities Data in R.

Arnold and Tilton 2015 Notes

Scatterplots: -scatter plot: thing of x= longitude, y= latitude To make scatter plot with underlying map: 1. plot your points 2. snippets package (maybe can also use OpenStreetMap package) pulls in an underlying map of ~25 tiles using osmap 3. then use points function to put the points back ontop of the map tiles 4. to set size of map, in plot, asp parameter sets the aspect ratio,the ratio of the scale on the y-axis to the scale of the x-axis, of a plot. If this is set to the ratio of the length of a degree of latitude to the length of a degree of longitude, the resulting map will be undistorted. -maps from Stamen Design with osmap function are black & white and less busy than OpenStreetMaps

ShapeFiles: -Map tiles don’t actually have information about the roads or geo features, they are just pixels (raster data.) -Vector data has meta-data attached to objects, and shapefiles are a way to store vectorize geospatial data -the SpatialPolygonsDataFrame has both geospatial information and a dataframe with metadata -plot these with projections, use rgdal package if needed to convert projections

Then read through ‘Spatial Data Manipulation’ on RSpatial. [http://rspatial.org/spatial/index.html]. There’s a lot here, so feel free to move through it quickly and think of this as reference.

1. From scatterplots to maps: Create a series of point maps with context

#a) Provided Map
datapath <- "../GEO200CN/data"
UCplaces <- read.csv(file.path(datapath, "UCplaces.csv"))
map <- openmap (c(38.55, -121.77), c(38.53, -121.74)) #openmap function gets a map using upper-left and lower-right corners (function knows that order)
plot (map) 

mapLatLong <- openproj(map) #openproj converts from the default mercator projection to Latitude Longitude 
plot (mapLatLong) 
points (UCplaces$long, UCplaces$lat, pch=16) #add the points 
text (UCplaces$long, UCplaces$lat, labels = UCplaces$names, pos=3) #add  labels

# b) change the background context map, [hint: ?openmap will give you a list of possible types]
map2 <- openmap(c(38.55, -121.77), c(38.53, -121.74), zoom = NULL, type = "stamen-watercolor")
plot(map2)

map2LatLong <- openproj(map2)
plot(map2LatLong)
points (UCplaces$long, UCplaces$lat, pch=16) 
text (UCplaces$long, UCplaces$lat, labels = UCplaces$names, pos=3)

#c) change the extent of the map to include more of Davis, and change the size and appearance of the points and labels as appropriate
map3 <- openmap(c(38.562083, -121.787334), c(38.526302, -121.728406), zoom = NULL, type = "esri")
plot(map3)

map3LatLong <- openproj(map3)
plot(map3LatLong)
points (UCplaces$long, UCplaces$lat, col= "blue", bg="blue", cex= 1.5, pch=23) 
text(UCplaces$long, UCplaces$lat, labels = UCplaces$names, col = "blue", cex = 1.2, pos=4)

2. Map displays: Create a series of simple maps with data

#a) Lat Long map of all the states and Puerto Rico
#read in the shape file
datapath1 <- "../GEO200CN/data/State_2010Census_DP1"
statefp <- file.path(datapath1,"State_2010Census_DP1")
states <-readOGR(dsn = statefp, layer = "State_2010Census_DP1" ) 
## OGR data source with driver: ESRI Shapefile 
## Source: "../GEO200CN/data/State_2010Census_DP1/State_2010Census_DP1", layer: "State_2010Census_DP1"
## with 52 features
## It has 195 fields
#readOGR reads a shapefile in, need to define first argument, DSN, is the path or director, seond argument, layer, is the file itself

class(states)
## [1] "SpatialPolygonsDataFrame"
## attr(,"package")
## [1] "sp"
dim(states) #52 rows which include Washington DC, and PR
## [1]  52 195
index <- (as.data.frame(states)$STUSPS10 %in% c("AK", "HI")) #removes AK and HI from the map
states1 <- states[!index,]
dim(states1)
## [1]  50 195
class(states1)
## [1] "SpatialPolygonsDataFrame"
## attr(,"package")
## [1] "sp"
plot(states1)

#identify the projection this is in, then it can be converted if necessary
crs(states1) # +proj=longlat +datum=NAD83 +no_defs +ellps=GRS80 +towgs84=0,0,0 
## CRS arguments:
##  +proj=longlat +datum=NAD83 +no_defs +ellps=GRS80 +towgs84=0,0,0
# b) a map in the UTM coordinate system with the states labeled
#install.packages("rgeos")
library(rgeos)
## rgeos version: 0.3-23, (SVN revision 546)
##  GEOS runtime version: 3.5.0-CAPI-1.9.0 r4084 
##  Linking to sp version: 1.2-4 
##  Polygon checking: TRUE
#need to assign a projection, then it can be transformed
#latlong<- CRS(projargs = "+proj=longlat")
#proj4string(states1) <- latlong
#Actually, this wasn't necessary, as latlong had already been assigned (which I knew b/c I checked above). Good to know how to assign this if necessary, here it wasn't I just need to get to the converting part

utms <- CRS(projargs = "+proj=utm +zone=14")
states1Trans <- spTransform(x = states1, CRSobj = utms) #transformed to utms
plot(states1Trans)

centroid <- gCentroid(spgeom=states1Trans, byid=TRUE) #find the middle of the states, gCentroid returns the center of the states
head(centroid)
## class       : SpatialPoints 
## features    : 1 
## extent      : -198256.2, -198256.2, 4796062, 4796062  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=utm +zone=14 +ellps=WGS84
text(x=centroid$x, 
     y=centroid$y,
     label=as.data.frame(states1Trans)$NAME10, cex=0.7)

#c) a cross hatch map showing the proportion of seasonal of recreational housing in each state

states1Transdata <-as.data.frame(states1)
#data now in data frame, now we look at the housing and how it is binned (seasonal, rec, or occaisional)

perHouseRec <- states1Transdata$DP0180008 / states1Transdata$DP0180001  #finds fancy housing/overall housing
bins <- quantile(perHouseRec, seq(0,1,length.out=8))  #quantile makes sample quantiles (weighted averages) to correspond to given probabilities. perHouseRec is the vector being sampled (all the fancy housing), and seq refers to what probabilities are being used. this makes 8 bins between 0 and 1

x <- 2:18
v <- c(5, 10, 15) # create two bins [5,10) and [10,15). actually, this makes four bins,[0,5], [5,10], [10,15], [15,]
findInterval(x, v) #tells how many of x fit into each of the 4 bins
##  [1] 0 0 0 1 1 1 1 1 2 2 2 2 2 3 3 3 3
binId <- findInterval(perHouseRec, bins) #findIntervalmeans find each interval of bins that contains perHouseRec
binId #this shows which bin each of the elements in perHouseRec should fall in, use this for below:
##  [1] 6 3 1 5 3 5 4 4 7 4 3 2 4 3 1 1 6 3 2 5 6 6 5 7 7 4 2 3 2 5 7 6 4 4 6
## [36] 7 5 1 2 1 1 7 6 2 1 3 2 7 8 5
densityVals <- seq_len(length(bins)) * 5  #seq_len generates regula sequences in the data bins, and is than scaled to a factor of 5
plot(states1Trans, density=densityVals[binId]) #plotting the raser object states1Trans, and a second object wich is a function of some data

#d) produce two additional maps using the color palettes defined above as colSet01 and colSet02
colSet01 = c("grey90", "grey80", "grey70", "grey60", "grey50", "grey40", "grey30", "grey20") 
plot (states1Trans, col=colSet01[binId]) #that's cool! uses colors to plot it, rather than the density hashmarks

colSet02 = rev(heat.colors(length (bins))) #rev reverses the argument, heat.colors is a series of colors
plot (states1Trans, col=colSet02[binId])

3. Exploring spatial data: Get to know the structure of spatial data in R

#a) Read the shapefile "galap.shp" 
galap.a <- file.path(datapath,"galap")
galap <-readOGR(dsn = galap.a, layer = "galap" )
## OGR data source with driver: ESRI Shapefile 
## Source: "../GEO200CN/data/galap", layer: "galap"
## with 30 features
## It has 12 fields
#b) Show the class of the object in R
classgalap <- class(galap)

#c) How many polygons are there in this object?
galap
## class       : SpatialPolygonsDataFrame 
## features    : 30 
## extent      : 610242.7, 918583.1, -156287.7, 185890.3  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=utm +zone=15 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0 
## variables   : 12
## names       :  NAME_0,     NAME_1,     NAME_2, Island, ID, species, native,    area, elevation, distNear, distSC, areaAdj 
## min values  : Ecuador, Galápagos,    Isabela, Baltra,  1,       2,      0,    0.01,         0,      0.2,    0.0,    0.03 
## max values  : Ecuador, Galápagos, Santa Cruz,   Wolf, 30,     444,     95, 4669.32,      1707,     47.4,  290.2, 4669.32
polygonsgalap <- 30

#d) How many variables are there in this object?
variablesgalap <- 12
  1. classes of galap = SpatialPolygonsDataFrame
  2. polygons in galap = 30
  3. variables = 12
#e) Make a plot of the islands (their outlines).
plot(galap)

#f) Make a scatter plot of the number of species as a function of the size of the island. 
galapdf <- as.data.frame(galap)
plot(x = galapdf$area, y = galapdf$species, pch = 21, bg = "blue")
m <- lm(species ~ area, data = galapdf) 
m
## 
## Call:
## lm(formula = species ~ area, data = galapdf)
## 
## Coefficients:
## (Intercept)         area  
##    63.78328      0.08196
abline(a= 63.78328, b= 0.08196)

  1. What quantity would you use to make a choropleth to represent the number of species on each island? A: I need a value of species/area for each island, then the colors will be attached to this value
#h) Create that quantity as a new variable in the SpatialPolyonsDataFrame and plot it with spplot 

galapdf <- transform(galapdf, sparea = species / area)
class(galapdf)
## [1] "data.frame"
galap1 <- SpatialPolygonsDataFrame(galap, data= galapdf, match.ID = TRUE)
#b/c galapdf was just a data.frame, I had to convert it back to a spatial object in order to use SpatialPolygonsDataFrame. for this function, the first argument is that object with the polygons (galap) and the second argument is the data frame with attributes that I want combined with it

spplot(galap1, "sparea", main = "Galapagos Species Per Island Area" )

#for spplot, 1st argument is the data, 2nd argument in quotes is the variable I want to highlight, and I added main, which is the title

#i) Select the largest island and save it to a new shapefile.

isabelasp <- galap1[2,]
#[] these brackets subset spatial (or other) data. so, here, we are just pulling out the row that has the info for isabela island
plot(isabelasp)

class(isabelasp)
## [1] "SpatialPolygonsDataFrame"
## attr(,"package")
## [1] "sp"
writeOGR(isabelasp, dsn = datapath, layer = "isabelasp", driver = "ESRI Shapefile" ) #this works just like readOGR, 1st argument is the data, dsn is your file path, layer is the name of the shapefile. the new thing is adding a driver, which says how you want the shapefile out
# j) Download elevation data for Ecuador (use the function getData in the raster package)
getData('ISO3')
##     ISO3                                         NAME
## 1    ABW                                        Aruba
## 2    AFG                                  Afghanistan
## 3    AGO                                       Angola
## 4    AIA                                     Anguilla
## 5    ALA                                        Åland
## 6    ALB                                      Albania
## 7    AND                                      Andorra
## 8    ARE                         United Arab Emirates
## 9    ARG                                    Argentina
## 10   ARM                                      Armenia
## 11   ASM                               American Samoa
## 12   ATA                                   Antarctica
## 13   ATF                  French Southern Territories
## 14   ATG                          Antigua and Barbuda
## 15   AUS                                    Australia
## 16   AUT                                      Austria
## 17   AZE                                   Azerbaijan
## 18   BDI                                      Burundi
## 19   BEL                                      Belgium
## 20   BEN                                        Benin
## 21   BES            Bonaire, Saint Eustatius and Saba
## 22   BFA                                 Burkina Faso
## 23   BGD                                   Bangladesh
## 24   BGR                                     Bulgaria
## 25   BHR                                      Bahrain
## 26   BHS                                      Bahamas
## 27   BIH                       Bosnia and Herzegovina
## 28   BLM                             Saint-Barthélemy
## 29   BLR                                      Belarus
## 30   BLZ                                       Belize
## 31   BMU                                      Bermuda
## 32   BOL                                      Bolivia
## 33   BRA                                       Brazil
## 34   BRB                                     Barbados
## 35   BRN                                       Brunei
## 36   BTN                                       Bhutan
## 37   BVT                                Bouvet Island
## 38   BWA                                     Botswana
## 39   CAF                     Central African Republic
## 40   CAN                                       Canada
## 41   CCK                                Cocos Islands
## 42   CHE                                  Switzerland
## 43   CHL                                        Chile
## 44   CHN                                        China
## 45   CIV                                Côte d'Ivoire
## 46   CMR                                     Cameroon
## 47   COD             Democratic Republic of the Congo
## 48   COG                            Republic of Congo
## 49   COK                                 Cook Islands
## 50   COL                                     Colombia
## 51   COM                                      Comoros
## 52   CPV                                   Cape Verde
## 53   CRI                                   Costa Rica
## 54   CUB                                         Cuba
## 55   CUW                                      Curaçao
## 56   CXR                             Christmas Island
## 57   CYM                               Cayman Islands
## 58   CYP                                       Cyprus
## 59   CZE                               Czech Republic
## 60   DEU                                      Germany
## 61   DJI                                     Djibouti
## 62   DMA                                     Dominica
## 63   DNK                                      Denmark
## 64   DOM                           Dominican Republic
## 65   DZA                                      Algeria
## 66   ECU                                      Ecuador
## 67   EGY                                        Egypt
## 68   ERI                                      Eritrea
## 69   ESH                               Western Sahara
## 70   ESP                                        Spain
## 71   EST                                      Estonia
## 72   ETH                                     Ethiopia
## 73   FIN                                      Finland
## 74   FJI                                         Fiji
## 75   FLK                             Falkland Islands
## 76   FRA                                       France
## 77   FRO                                Faroe Islands
## 78   FSM                                   Micronesia
## 79   GAB                                        Gabon
## 80   GBR                               United Kingdom
## 81   GEO                                      Georgia
## 82   GGY                                     Guernsey
## 83   GHA                                        Ghana
## 84   GIB                                    Gibraltar
## 85   GIN                                       Guinea
## 86   GLP                                   Guadeloupe
## 87   GMB                                       Gambia
## 88   GNB                                Guinea-Bissau
## 89   GNQ                            Equatorial Guinea
## 90   GRC                                       Greece
## 91   GRD                                      Grenada
## 92   GRL                                    Greenland
## 93   GTM                                    Guatemala
## 94   GUF                                French Guiana
## 95   GUM                                         Guam
## 96   GUY                                       Guyana
## 97   HKG                                    Hong Kong
## 98   HMD            Heard Island and McDonald Islands
## 99   HND                                     Honduras
## 100  HRV                                      Croatia
## 101  HTI                                        Haiti
## 102  HUN                                      Hungary
## 103  IDN                                    Indonesia
## 104  IMN                                  Isle of Man
## 105  IND                                        India
## 106  IOT               British Indian Ocean Territory
## 107  IRL                                      Ireland
## 108  IRN                                         Iran
## 109  IRQ                                         Iraq
## 110  ISL                                      Iceland
## 111  ISR                                       Israel
## 112  ITA                                        Italy
## 113  JAM                                      Jamaica
## 114  JEY                                       Jersey
## 115  JOR                                       Jordan
## 116  JPN                                        Japan
## 117  KAZ                                   Kazakhstan
## 118  KEN                                        Kenya
## 119  KGZ                                   Kyrgyzstan
## 120  KHM                                     Cambodia
## 121  KIR                                     Kiribati
## 122  KNA                        Saint Kitts and Nevis
## 123  KOR                                  South Korea
## 124  KWT                                       Kuwait
## 125  LAO                                         Laos
## 126  LBN                                      Lebanon
## 127  LBR                                      Liberia
## 128  LBY                                        Libya
## 129  LCA                                  Saint Lucia
## 130  LIE                                Liechtenstein
## 131  LKA                                    Sri Lanka
## 132  LSO                                      Lesotho
## 133  LTU                                    Lithuania
## 134  LUX                                   Luxembourg
## 135  LVA                                       Latvia
## 136  MAC                                        Macao
## 137  MAF                                 Saint-Martin
## 138  MAR                                      Morocco
## 139  MCO                                       Monaco
## 140  MDA                                      Moldova
## 141  MDG                                   Madagascar
## 142  MDV                                     Maldives
## 143  MEX                                       Mexico
## 144  MHL                             Marshall Islands
## 145  MKD                                    Macedonia
## 146  MLI                                         Mali
## 147  MLT                                        Malta
## 148  MMR                                      Myanmar
## 149  MNE                                   Montenegro
## 150  MNG                                     Mongolia
## 151  MNP                     Northern Mariana Islands
## 152  MOZ                                   Mozambique
## 153  MRT                                   Mauritania
## 154  MSR                                   Montserrat
## 155  MTQ                                   Martinique
## 156  MUS                                    Mauritius
## 157  MWI                                       Malawi
## 158  MYS                                     Malaysia
## 159  MYT                                      Mayotte
## 160  NAM                                      Namibia
## 161  NCL                                New Caledonia
## 162  NER                                        Niger
## 163  NFK                               Norfolk Island
## 164  NGA                                      Nigeria
## 165  NIC                                    Nicaragua
## 166  NIU                                         Niue
## 167  NLD                                  Netherlands
## 168  NOR                                       Norway
## 169  NPL                                        Nepal
## 170  NRU                                        Nauru
## 171  NZL                                  New Zealand
## 172  OMN                                         Oman
## 173  PAK                                     Pakistan
## 174  PAN                                       Panama
## 175  PCN                             Pitcairn Islands
## 176  PER                                         Peru
## 177  PHL                                  Philippines
## 178  PIS                              Paracel Islands
## 179  PLW                                        Palau
## 180  PNG                             Papua New Guinea
## 181  POL                                       Poland
## 182  PRI                                  Puerto Rico
## 183  PRK                                  North Korea
## 184  PRT                                     Portugal
## 185  PRY                                     Paraguay
## 186  PSE                                    Palestina
## 187  PYF                             French Polynesia
## 188  QAT                                        Qatar
## 189  REU                                      Reunion
## 190  ROU                                      Romania
## 191  RUS                                       Russia
## 192  RWA                                       Rwanda
## 193  SAU                                 Saudi Arabia
## 194  SDN                                        Sudan
## 195  SEN                                      Senegal
## 196  SGP                                    Singapore
## 197  SGS South Georgia and the South Sandwich Islands
## 198  SHN                                 Saint Helena
## 199  SJM                       Svalbard and Jan Mayen
## 200  SLB                              Solomon Islands
## 201  SLE                                 Sierra Leone
## 202  SLV                                  El Salvador
## 203  SMR                                   San Marino
## 204  SOM                                      Somalia
## 205  SP-                              Spratly islands
## 206  SPM                    Saint Pierre and Miquelon
## 207  SRB                                       Serbia
## 208  SSD                                  South Sudan
## 209  STP                        Sao Tome and Principe
## 210  SUR                                     Suriname
## 211  SVK                                     Slovakia
## 212  SVN                                     Slovenia
## 213  SWE                                       Sweden
## 214  SWZ                                    Swaziland
## 215  SXM                                 Sint Maarten
## 216  SYC                                   Seychelles
## 217  SYR                                        Syria
## 218  TCA                     Turks and Caicos Islands
## 219  TCD                                         Chad
## 220  TGO                                         Togo
## 221  THA                                     Thailand
## 222  TJK                                   Tajikistan
## 223  TKL                                      Tokelau
## 224  TKM                                 Turkmenistan
## 225  TLS                                   East Timor
## 226  TON                                        Tonga
## 227  TTO                          Trinidad and Tobago
## 228  TUN                                      Tunisia
## 229  TUR                                       Turkey
## 230  TUV                                       Tuvalu
## 231  TWN                                       Taiwan
## 232  TZA                                     Tanzania
## 233  UGA                                       Uganda
## 234  UKR                                      Ukraine
## 235  UMI         United States Minor Outlying Islands
## 236  URY                                      Uruguay
## 237  USA                                United States
## 238  UZB                                   Uzbekistan
## 239  VAT                                 Vatican City
## 240  VCT             Saint Vincent and the Grenadines
## 241  VEN                                    Venezuela
## 242  VGB                       British Virgin Islands
## 243  VIR                         Virgin Islands, U.S.
## 244  VNM                                      Vietnam
## 245  VUT                                      Vanuatu
## 246  WLF                            Wallis and Futuna
## 247  WSM                                        Samoa
## 248  XAD                        Akrotiri and Dhekelia
## 249  XCA                                  Caspian Sea
## 250  XCL                            Clipperton Island
## 251  XKO                                       Kosovo
## 252  XNC                              Northern Cyprus
## 253  YEM                                        Yemen
## 254  ZAF                                 South Africa
## 255  ZMB                                       Zambia
## 256  ZWE                                     Zimbabwe
ecualt <-getData("alt", country = "ECU", mask = TRUE)

# k) Use the crop function in the raster package and then map the elevation data (add the island outlines).

#When I tried to crop these together, they didn't overlap in extent, maybe b/c they are projecting in different systems?
crs(ecualt)
## CRS arguments:
##  +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0
crs(galap)
## CRS arguments:
##  +proj=utm +zone=15 +datum=WGS84 +units=m +no_defs +ellps=WGS84
## +towgs84=0,0,0
#transform altitude to utms. have to do this with galap bc ecualt is a raster object
long_lat <- CRS(projargs = " +proj=longlat +datum=WGS84")
galap_trans <- spTransform(x = galap, CRSobj = long_lat)

galapalt <- crop(ecualt, galap_trans)

#plot(galapalt) only plots the altitude
#plot(galap, add = TRUE) didn't add this

#From: http://rspatial.org/spatial/rst/9-maps.html. when the above doesn't work, use spplot. here in bounds, I pulled out the polygons from galap to be the boundaries, then plotted them wih the altitude
bounds <- list("sp.polygons", galap)
spplot(galapalt, sp.layout=bounds )

# l) Create a contour map of elevation in Ecuador.
galap_cont <- contour(galapalt) #that's fine, but it's not amazing. it looks like ther is a lot of other cool stuff to do with contours, too.

#spplot(galapalt, sp.layout=galap_cont)
#plot(galapalt)
#plot(galap_cont, add = TRUE)

#another good mapping page: http://www.nickeubank.com/wp-content/uploads/2015/10/RGIS3_MakingMaps_part1_mappingVectorData.html#spplot


#contour(galapalt, add = TRUE) this didn't work, I think different projection systems?

filledContour(galapalt)

#another contour map option

Upload to smartsite an Rmarkdown and HTML file.

Due: Monday, April 24, 9a